library(survey)
library(dplyr)
library(stargazer)
library(performance)
library(jtools)
library(ggplot2)
library(gridExtra)

# Code variables of interest as binaries for wave 1
summary(factor(mil_conf.df$Q15A))

mil_conf.df$Q15A2[mil_conf.df$Q15A < 77] <- 0
mil_conf.df$Q15A2[mil_conf.df$Q15A < 3] <- 1
mil_conf.df$Q15B2[mil_conf.df$Q15B < 77] <- 0
mil_conf.df$Q15B2[mil_conf.df$Q15B < 3] <- 1
mil_conf.df$Q15C2[mil_conf.df$Q15C < 77] <- 0
mil_conf.df$Q15C2[mil_conf.df$Q15C < 3] <- 1
mil_conf.df$Q15D2[mil_conf.df$Q15D < 77] <- 0
mil_conf.df$Q15D2[mil_conf.df$Q15D < 3] <- 1
mil_conf.df$Q15E2[mil_conf.df$Q15E < 77] <- 0
mil_conf.df$Q15E2[mil_conf.df$Q15E < 3] <- 1
mil_conf.df$Q15F2[mil_conf.df$Q15F < 77] <- 0
mil_conf.df$Q15F2[mil_conf.df$Q15F < 3] <- 1
mil_conf.df$Q15G2[mil_conf.df$Q15G < 77] <- 0
mil_conf.df$Q15G2[mil_conf.df$Q15G < 3] <- 1

mil_conf.df$Q13T <- 0
mil_conf.df$Q13T[mil_conf.df$Q13 < 3] <- 1

# Code variables of interest as binaries for wave 2
# Institutions
df$Q8A2[df$Q8A < 77] <- 0
df$Q8A2[df$Q8A < 3] <- 1
df$Q8B2[df$Q8B < 77] <- 0
df$Q8B2[df$Q8B < 3] <- 1
df$Q8C2[df$Q8C < 77] <- 0
df$Q8C2[df$Q8C < 3] <- 1
df$Q8D2[df$Q8D < 77] <- 0
df$Q8D2[df$Q8D < 3] <- 1
df$Q8E2[df$Q8E < 77] <- 0
df$Q8E2[df$Q8E < 3] <- 1
df$Q8F2[df$Q8F < 77] <- 0
df$Q8F2[df$Q8F < 3] <- 1
df$Q8G2[df$Q8G < 77] <- 0
df$Q8G2[df$Q8G < 3] <- 1
df$Q8H2[df$Q8H < 77] <- 0
df$Q8H2[df$Q8H < 3] <- 1

summary(factor(df$Q7))
df$Q7T[df$Q7 < 77] <- 0
df$Q7T[df$Q7 < 3] <- 1
summary(factor(df$Q7T))

# Professions questions
df$Q10A2[df$Q10A < 77] <- 0
df$Q10A2[df$Q10A < 3] <- 1
df$Q10B2[df$Q10B < 77] <- 0
df$Q10B2[df$Q10B < 3] <- 1
df$Q10C2[df$Q10C < 77] <- 0
df$Q10C2[df$Q10C < 3] <- 1
df$Q10D2[df$Q10D < 77] <- 0
df$Q10D2[df$Q10D < 3] <- 1
df$Q10E2[df$Q10E < 77] <- 0
df$Q10E2[df$Q10E < 3] <- 1
df$Q10F2[df$Q10F < 77] <- 0
df$Q10F2[df$Q10F < 3] <- 1
df$Q10G2[df$Q10G < 77] <- 0
df$Q10G2[df$Q10G < 3] <- 1
df$Q10H2[df$Q10H < 77] <- 0
df$Q10H2[df$Q10H < 3] <- 1

summary(factor(df$Q9))
df$Q9T[df$Q9 < 77] <- 0
df$Q9T[df$Q9 < 3] <- 1
summary(factor(df$Q9T))



# Creating weighted survey design objects
w1_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

w2_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight2,
    data = df
  )

# Creating Figure 4.9
# Wave 1
m1a <- svyglm(Q13T ~ Q15A2 + Q15B2 + Q15C2 + Q15D2 + Q15E2 + Q15F2 + Q15G2 + 
              Q18a + dem + rep + ideo3 + male + white + black + hispanic + asian +
              income5 + educ5 + unemployed + married + silent + boomer + genx + millen +
              activeduty + vet + social + family + midwest + south + west + catholic + christian + norelig +
              city + rural,
             family = binomial(link = "logit"),
             design = w1_design)
summary(m1a)

m1b <- svyglm(Q13T ~ Q15A2 + Q15B2 + Q15C2 + Q15D2 + Q15E2 + Q15F2 + Q15G2 + 
                dem + rep + ideo3 + male + white + black + hispanic + asian +
                income5 + educ5 + unemployed + married + silent + boomer + genx + millen +
                activeduty + vet + social + family + midwest + south + west + catholic + christian + norelig +
                city + rural,
              family = binomial(link = "logit"),
              design = w1_design)
summary(m1b)

# Wave 2
# Jim creates INSTITUTION variable and PROFESSION variable.

m2 <- svyglm(Q7T ~ Q8A2 + Q8B2 + Q8C2 + Q8D2 + Q8E2 + Q8F2 + Q8G2 + 
               dem + rep + ideo3 + male + white + black + hispanic + asian +
               income5 + EDUC5 + unemployed + married + silent + boomer + genx + millen +
               activeduty + vet + social + family + midwest + south + west + catholic + christian + norelig +
               city + rural,
             family = binomial(link = "logit"),
             design = w2_design)

m3 <- svyglm(Q9T ~ Q10A2 + Q10B2 + Q10C2 + Q10D2 + Q10E2 + Q10F2 + Q10G2 + 
               dem + rep + ideo3 + male + white + black + hispanic + asian +
               income5 + EDUC5 + unemployed + married + silent + boomer + genx + millen +
               activeduty + vet + social + family + midwest + south + west + catholic + christian + norelig +
               city + rural,
             family = binomial(link = "logit"),
             design = w2_design)


p1 <- plot_summs(m1b, legend.pos = "none") + 
  xlim(-2,2) + labs(title="Institutions (2019)") +
  theme(axis.title.x = element_blank(), plot.title = element_text(hjust=0.5))
p2 <- plot_summs(m2, legend.pos = "none") + 
  xlim(-2,2) + labs(title="Agencies (2020)") + 
  theme(axis.title.x = element_blank(), plot.title = element_text(hjust=0.5))
p3 <- plot_summs(m3, legend.pos = "none") + 
  xlim(-2,2) + labs(title="Professions (2020)") +
  theme(axis.title.x = element_blank(), plot.title = element_text(hjust=0.5))

grid.arrange(p1,p2,p3,ncol=3)
